Question 1: How is grazer density related to macrohabitat composition?
Question 2: How are vegetation structure and functional traits influenced by grazing?
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
## Loading required package: spam
## Spam version 2.10-0 (2023-10-23) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: viridisLite
##
## Try help(fields) to get started.
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following object is masked from 'package:spam':
##
## det
## Loading required package: lme4
##
## arm (Version 1.14-4, built: 2024-4-1)
## Working directory is /Users/lillalovasz/Library/CloudStorage/Dropbox/papers/vegetation-grazer/rcode
##
## Attaching package: 'arm'
## The following object is masked from 'package:spam':
##
## display
#read data
dhorse19_20 <- read.table("../data/horses-all-positions_2019-jan-09_2020_aug-31.txt", header=TRUE, sep = "\t")
dhorse20_21 <- read.table("../data/horses-all-positions_2020-sept-01_2021_aug_31.txt", header=TRUE, sep = "\t")
dcattle <- read.table("../data/cattle_all-positions-2019-jan-09_2021-aug-31.txt", header=TRUE, sep = "\t")
dfencefull <- read.csv("../data/fence.csv", header=TRUE)
dfence <- read.table("../data/fence-reduced.txt", header=TRUE, sep="\t")
dfence <- dfencefull # when we include the forest
dcells <- read.table("../data/cells_info_veg-plots-with-forest.txt", header=TRUE, sep="\t")
dlandcover <- read.table("../data/habitat-with-forest.txt",header=TRUE, sep="\t")
dveg <- read.table("../data/veg_data_2018-2019-2020-2021-summer.txt",header=TRUE, sep="\t")Make all coordinates being the same
Clip grazer positions to the fence-area including the forest
Correct the number of locations for the duration of the season and nr of animals tracked.
## Name ID Date Time X Latitude
## 1 Onissoi 401776 1/9/2019 00:00 NA 47.62875
## 2 Odissa 401777 1/9/2019 00:00 NA 47.62950
## 3 TERMINATED - 401775 401775 1/9/2019 00:00 NA 47.62876
## 4 Odaline 401778 1/9/2019 00:00 NA 47.62933
## 5 TERMINATED - 401774 lost (ex-Olocen) 401774 1/9/2019 00:00 NA 47.62883
## 6 Odaline 401778 1/9/2019 01:00 NA 47.62976
## Longitude Speed
## 1 7.556987 0.00
## 2 7.556098 0.00
## 3 7.556587 0.37
## 4 7.556327 0.00
## 5 7.556852 2.04
## 6 7.556088 0.74
## 'data.frame': 102216 obs. of 8 variables:
## $ Name : chr "Onissoi" "Odissa" "TERMINATED - 401775" "Odaline" ...
## $ ID : int 401776 401777 401775 401778 401774 401778 401776 401774 401777 401775 ...
## $ Date : chr "1/9/2019" "1/9/2019" "1/9/2019" "1/9/2019" ...
## $ Time : chr "00:00" "00:00" "00:00" "00:00" ...
## $ X : logi NA NA NA NA NA NA ...
## $ Latitude : num 47.6 47.6 47.6 47.6 47.6 ...
## $ Longitude: num 7.56 7.56 7.56 7.56 7.56 ...
## $ Speed : num 0 0 0.37 0 2.04 0.74 0.19 2.22 0.56 0.19 ...
## [1] ""
## [2] "Odaline"
## [3] "Odissa"
## [4] "Onissoi"
## [5] "TERMINATED - 401707"
## [6] "TERMINATED - 401774 lost (ex-Olocen)"
## [7] "TERMINATED - 401775"
## [1] 27244
## Name ID Date Time X Latitude Longitude X.1 Speed
## 1 TERMINATED - 401773 401773 1/9/2019 06:00 NA 47.62071 7.539350 NA 0.74
## 2 TERMINATED - 401771 401771 1/9/2019 06:01 NA 47.62059 7.538698 NA 0.19
## [1] "401767" "401769" "401770" "401771" "401773" "421430" "421431"
## [8] "T5H-6839" "T5H-6840"
## [1] " TERMINATED - 401773" "Gevaudan (6759) Tellus" "Kayla (3412) Tellus"
## [4] "TERMINATED - 401767" "TERMINATED - 401769" "TERMINATED - 401770"
## [7] "TERMINATED - 401771" "TERMINATED - 421430" "TERMINATED - 421431"
dcattle$X <- NULL
dcattle$X.1 <- NULL
##
#merge horse and cattle
dgrazer <- rbind(dhorse, dcattle)
#create "Animal" type
dgrazer <- mutate(dgrazer, Animal = ifelse(dgrazer$ID %in% c("401707","401774","401775","401776","401777","401778"), "Horse", "Cattle"))
head(dgrazer,2)## Name ID Date Time Latitude Longitude Speed Animal
## 1 Onissoi 401776 1/9/2019 00:00 47.62875 7.556987 0 Horse
## 2 Odissa 401777 1/9/2019 00:00 47.62950 7.556098 0 Horse
## [1] 0
dgrazer$Animal <- as.factor(dgrazer$Animal)
#date in correct format
index <- grepl("/", dgrazer$Date)
dgrazer$Date[index] <- as.character(as.Date(as.character(dgrazer$Date[index]), format="%m/%d/%Y"))
dgrazer$Date[!index] <- as.character(as.Date(as.character(dgrazer$Date[!index]), format="%d.%m.%y"))
#date for creating yday (day of year as decimal nr)
dgrazer$DateTime <- as.POSIXct(paste(as.character(dgrazer$Date), as.character(dgrazer$Time), sep=" "),format="%Y-%m-%d %H:%M", tz="CET")
##linear timeline as a number passed from 2019-jan-1.
dgrazer$dayofyear <- strptime(dgrazer$Date, format="%Y-%m-%d")$yday
dgrazer$year <- strptime(dgrazer$Date, format="%Y-%m-%d")$year+1900 #
head(dgrazer,2)## Name ID Date Time Latitude Longitude Speed Animal DateTime
## 1 Onissoi 401776 2019-01-09 00:00 47.62875 7.556987 0 Horse 2019-01-09
## 2 Odissa 401777 2019-01-09 00:00 47.62950 7.556098 0 Horse 2019-01-09
## dayofyear year
## 1 8 2019
## 2 8 2019
## [1] 366
#coordinates in correct format
tmp <- dgrazer
coordinates(tmp) <- c("Longitude", "Latitude") # make into a spatial point dataframe (WGS84 = long /y/; lat /x/)
crs.geo <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(tmp) <- crs.geo #define projection
tmp2 <- spTransform(tmp, CRS("+init=epsg:32232")) # utm transformation## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
#dgrazer$Latitude.UTM <- coordinates(tmp2)[,"Latitude"]
#dgrazer$Longitude.UTM <- coordinates(tmp2)[,"Longitude"]
#
dgrazer$Latitude.UTM <- coordinates(tmp2)[,"coords.x2"]
dgrazer$Longitude.UTM <- coordinates(tmp2)[,"coords.x1"]
tmp3 <- dfence
coordinates(tmp3) <- c("x", "y") # make into a spatial point dataframe
crs.geo <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(tmp3) <- crs.geo #define projection
tmp4 <- spTransform(tmp3, CRS("+init=epsg:32232")) # utm transformation
dfence$Latitude.UTM <- coordinates(tmp4)[,"coords.x2"]
dfence$Longitude.UTM <- coordinates(tmp4)[,"coords.x1"]
plot(dfence$x, dfence$y)
text(dfence$x, dfence$y, labels=dfence$Name, cex=0.6)## Name y x Latitude.UTM Longitude.UTM
## 1 Marker 1 47.62670 7.559454 5275811 391758.2
## 2 Marker 2 47.62683 7.559599 5275826 391769.4
## 3 Marker 3 47.62706 7.559288 5275852 391746.5
## 4 Marker 4 47.62732 7.558995 5275881 391725.1
## 5 Marker 5 47.62760 7.558683 5275913 391702.2
## 6 Marker 6 47.62788 7.558317 5275945 391675.3
########
#track and delete duplications
dgrazer$DateTime <- paste(as.character(dgrazer$Date), dgrazer$Time, dgrazer$Name)
table(table(as.character(dgrazer$DateTime)))##
## 1 2 3 4 5
## 147695 363 79 11 5
index2 <- duplicated(dgrazer$DateTime)
dgrazer <- dgrazer[!index2,]
#plot raw grazer points (outside the grid too)
plot(dcells$x, dcells$y, type="n", xlim=c(min(dcells$x), max(dcells$x)+50),
ylim=c(min(dcells$y), max(dcells$y)+50))
rect(dcells$x, dcells$y, dcells$x+50, dcells$y+50, border = "grey")
text(dcells$x, dcells$y, labels=dcells$cellnr, cex=0.6)
points(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, col=alpha(ifelse(dgrazer$Animal=='Horse', "blue","red"), 0.3), cex=0.6, pch=16)
points(dfence$Longitude.UTM, dfence$Latitude.UTM)#points(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, col="orange", cex=0.6, pch=16)
#only positions inside the fence
dgrazer$Inside <- point.in.polygon(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, dfence$Longitude.UTM, dfence$Latitude.UTM)
dgrazer <- dgrazer[dgrazer$Inside==1,]
#all-fence-area
#plot animal position in all the enclosure
plot(dfence$Longitude.UTM, dfence$Latitude.UTM)
points(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, col=alpha(ifelse(dgrazer$Animal=='Horse', "blue","red"), 0.3), cex=0.6, pch=16)# grazer positions inside the grid cells
allexistingx <- unique(dcells$x)
allexistingy <- unique(dcells$y)
for(i in 1:nrow(dgrazer)){
index <- (dgrazer$Longitude.UTM[i]-allexistingx) >0 & (dgrazer$Longitude.UTM[i]-allexistingx) <=50
if(sum(index)>0) dgrazer$x.cell[i] <- allexistingx[index]
index <- (dgrazer$Latitude.UTM[i]-allexistingy) >0 & (dgrazer$Latitude.UTM[i]-allexistingy) <=50
if(sum(index)>0) dgrazer$y.cell[i] <- allexistingy[index]
}
dgrazer <- dgrazer[!is.na(dgrazer$y.cell),]
#get the cellnr in the dataframe
dgrazer$cellnr <- dcells$cellnr[match(paste(dgrazer$x.cell, dgrazer$y.cell), paste(dcells$x, dcells$y))]
length(unique(dgrazer$cellnr)) ## [1] 160
#
dgrazer <- dgrazer[!is.na(dgrazer$cellnr),] # -> delete positions outside the grid
###replot
plot(dcells$x, dcells$y, type="n", xlim=c(min(dcells$x), max(dcells$x)+50),
ylim=c(min(dcells$y), max(dcells$y)+50))
rect(dcells$x, dcells$y, dcells$x+50, dcells$y+50, border = "grey")
text(dcells$x, dcells$y, labels=dcells$cellnr, cex=0.6)
points(dgrazer$Longitude.UTM, dgrazer$Latitude.UTM, cex=0.6,
pch=16, col=rgb(1,1,0,0.2))Question 1 (macrohabitat use of grazers): we use the average number of locations within each cell per season and year, averaged over the individuals tracked and averaged per day
dgrazer$Date.date <- strptime(dgrazer$Date, format="%Y-%m-%d")
dgrazerforveg <- dgrazer # make a copy for the micro-analyses later
dgrazer$season[dgrazer$Date.date >= "2019-06-24" & dgrazer$Date.date <= "2019-08-31"] <- "summer"
dgrazer$year.season[dgrazer$Date.date >= "2019-06-24" & dgrazer$Date.date <= "2019-08-31"] <- 2019
dgrazer$duration[dgrazer$Date.date >= "2019-06-24" & dgrazer$Date.date <= "2019-08-31"] <- 7+31+31
dgrazer$season[dgrazer$Date.date >= "2020-06-01" & dgrazer$Date.date <= "2020-08-31"] <- "summer"
dgrazer$year.season[dgrazer$Date.date >= "2020-06-01" & dgrazer$Date.date <= "2020-08-31"] <- 2020
dgrazer$duration[dgrazer$Date.date >= "2020-06-01" & dgrazer$Date.date <= "2020-08-31"] <- 30+31+31
dgrazer$season[dgrazer$Date.date >= "2021-06-28" & dgrazer$Date.date <= "2021-08-31"] <- "summer"
dgrazer$duration[dgrazer$Date.date >= "2021-06-28" & dgrazer$Date.date <= "2021-08-31"] <- 3+31+31
dgrazer$year.season[dgrazer$Date.date >= "2021-06-28" & dgrazer$Date.date <= "2021-08-31"] <- 2021
dgrazer$season[dgrazer$Date.date >= "2019-01-09" & dgrazer$Date.date <= "2019-02-28"] <- "winter"
dgrazer$duration[dgrazer$Date.date >= "2019-01-09" & dgrazer$Date.date <= "2019-02-28"] <- 23+28
dgrazer$year.season[dgrazer$Date.date >= "2019-01-09" & dgrazer$Date.date <= "2019-02-28"] <- 2019
dgrazer$season[dgrazer$Date.date >= "2019-12-01" & dgrazer$Date.date <= "2020-02-29"] <- "winter"
dgrazer$duration[dgrazer$Date.date >= "2019-12-01" & dgrazer$Date.date <= "2020-02-29"] <-
31+31+29
dgrazer$year.season[dgrazer$Date.date >= "2019-12-01" & dgrazer$Date.date <= "2020-02-29"] <- 2020
dgrazer$season[dgrazer$Date.date >= "2020-12-01" & dgrazer$Date.date <= "2021-02-28"] <- "winter"
dgrazer$duration[dgrazer$Date.date >= "2020-12-01" & dgrazer$Date.date <= "2021-02-28"] <- 31+31+28
dgrazer$year.season[dgrazer$Date.date >= "2020-12-01" & dgrazer$Date.date <= "2021-02-28"] <- 2021
dgrazer <- dgrazer[!is.na(dgrazer$season),]
dgrazer$year.season <- factor(dgrazer$year.season)
dgrazer$season <- factor(dgrazer$season)
dgrazer$cellnr <- factor(dgrazer$cellnr)
dgrazersum <- aggregate(Animal~year.season+season+cellnr, dgrazer, FUN=function(x) sum(x=="Cattle"), drop=FALSE)
names(dgrazersum)[names(dgrazersum)=="Animal"] <- "Cattle"
dgrazersum$Horse <- aggregate(Animal~year.season+season+cellnr, dgrazer, FUN=function(x) sum(x=="Horse"), drop=FALSE)$Animal
dgrazersum$Cattle[is.na(dgrazersum$Cattle)] <- 0
dgrazersum$Horse[is.na(dgrazersum$Horse)] <- 0
dgrazersum$duration <- dgrazer$duration[match(paste(dgrazersum$year.season, dgrazersum$season), paste(dgrazer$year.season, dgrazer$season))]
dgrazersum$grazer <- dgrazersum$Cattle + dgrazersum$Horse
dgrazersum$grazer <- dgrazersum$grazer/dgrazersum$duration
dgrazersum$Cattle <- dgrazersum$Cattle/dgrazersum$duration
dgrazersum$Horse <- dgrazersum$Horse/dgrazersum$duration
# number of individuals tracked
dtemp <- aggregate(ID~year+season, dgrazer[dgrazer$Animal=="Horse",], FUN=function(x) length(unique(x)), drop=FALSE)
dgrazersum$Horse_nrId <- dtemp$ID[match(paste(dgrazersum$year, dgrazersum$season),
paste(dtemp$year, dtemp$season))]
dtemp <- aggregate(ID~year+season, dgrazer[dgrazer$Animal=="Cattle",], FUN=function(x) length(unique(x)), drop=FALSE)
dgrazersum$Cattle_nrId <- dtemp$ID[match(paste(dgrazersum$year, dgrazersum$season),
paste(dtemp$year, dtemp$season))]
dgrazersum$horseUDcorr <- dgrazersum$Horse/dgrazersum$Horse_nrId
dgrazersum$cattleUDcorr <- dgrazersum$Cattle/dgrazersum$Cattle_nrId
dgrazersum$grazerUDcorr <- dgrazersum$horseUDcorr+dgrazersum$cattleUDcorr
hist(dgrazersum$grazerUDcorr)## [1] 0.000000 3.813043
dgrazersum$x <- dcells$x[match(dgrazersum$cellnr, dcells$cellnr)]
dgrazersum$y <- dcells$y[match(dgrazersum$cellnr, dcells$cellnr)]
dgrazersum$grazer.log <- log(dgrazersum$grazerUDcorr+min(dgrazersum$grazer[dgrazersum$grazer>0])/2)
dgrazersum$horse.log <- log(dgrazersum$horseUDcorr++min(dgrazersum$Horse[dgrazersum$Horse>0])/2)
dgrazersum$cattle.log <- log(dgrazersum$cattleUDcorr++min(dgrazersum$Cattle[dgrazersum$Cattle>0])/2)
maxloggrazer <- max(dgrazersum$grazer.log)
minloggrazer <- min(dgrazersum$grazer.log)
dgrazersum$colnum <- (dgrazersum$grazer.log-minloggrazer)/(maxloggrazer-minloggrazer)
par(mfrow=c(3,2), mar=c(1,1,2,1))
for(i in levels(dgrazersum$year.season)){
index <- dgrazersum$season=="winter"&dgrazersum$year.season==i
plot(dgrazersum$x[index], dgrazersum$y[index], pch=15, col=grey(1-dgrazersum$colnum[index]), cex=3, main=paste("Winter",i))
index <- dgrazersum$season=="summer"&dgrazersum$year.season==i
plot(dgrazersum$x[index], dgrazersum$y[index], pch=15, col=grey(1-dgrazersum$colnum[index]), cex=3, main=paste("Summer", i))
}##put landcover in the dataframe
#head(dlandcover)
dat <- merge(dgrazersum, dlandcover, by="cellnr")
#head(dat,2)
cov.varnames <- c("Cover.tree", "Cover.sapling", "Cover.shrub", "Cover.meadow", "Surface.water", "Bareground")
pairs(dat[, cov.varnames])Figure 2.1: Pairs-plot of the cover variables
Question 2 (how do vegetation traits change as a function of grazer density): we use the grazer density as a predictor for the vegetation structure. Therefore, we need to measure the total grazing density summed over all individuals. Because not all individuals were tracked, and because both grazer species are known to stick together in group, we multiplied the average locations per individual by the number of animals that were actually grazing.
# need to add in dveg the grazing density
dgrazerforveg$include <- FALSE
dgrazerforveg$include[dgrazerforveg$Date.date >= "2019-01-09" & dgrazerforveg$Date.date <= "2019-06-01"] <- TRUE
dgrazerforveg$year.season[dgrazerforveg$Date.date >= "2019-01-09" & dgrazerforveg$Date.date <= "2019-06-01"] <- 2019
dgrazerforveg$include[dgrazerforveg$Date.date >= "2020-01-01" & dgrazerforveg$Date.date <= "2020-06-01"] <- TRUE
dgrazerforveg$year.season[dgrazerforveg$Date.date >= "2020-01-01" & dgrazerforveg$Date.date <= "2020-06-01"] <- 2020
dgrazerforveg$include[dgrazerforveg$Date.date >= "2021-01-01" & dgrazerforveg$Date.date <= "2021-06-01"] <- TRUE
dgrazerforveg$year.season[dgrazerforveg$Date.date >= "2021-01-01" & dgrazerforveg$Date.date <= "2021-06-01"] <- 2021
dgrazerforveg <- dgrazerforveg[dgrazerforveg$include,]
dgrazersumveg <- aggregate(Animal~year.season+cellnr, dgrazerforveg, FUN=function(x) sum(x=="Cattle"), drop=FALSE)
names(dgrazersumveg)[names(dgrazersumveg)=="Animal"] <- "Cattle"
dgrazersumveg$Horse <- aggregate(Animal~year.season+cellnr, dgrazerforveg, FUN=function(x) sum(x=="Horse"), drop=FALSE)$Animal
dgrazersumveg$Cattle[is.na(dgrazersumveg$Cattle)] <- 0
dgrazersumveg$Horse[is.na(dgrazersumveg$Horse)] <- 0
# number of individuals tracked measured in "animal-days"
dtemp<- aggregate(dayofyear~Name + year, data=dgrazerforveg, function(x) length(unique(x)))
dtemp$Animal <- dgrazer$Animal[match(dtemp$Name, dgrazer$Name)]
dtemp1 <- aggregate(dayofyear~year+Animal, data=dtemp, sum)
dgrazersumveg$nrCattletracked <- dtemp1$dayofyear[dtemp1$Animal=="Cattle"][match(dgrazersumveg$year.season, dtemp1$year[dtemp1$Animal=="Cattle"])]
dgrazersumveg$nrHorsetracked <- dtemp1$dayofyear[dtemp1$Animal=="Horse"][match(dgrazersumveg$year.season, dtemp1$year[dtemp1$Animal=="Horse"])]
drealpres <- read.table("../data/animal-real-presence.txt", header=TRUE, sep="\t")
drealpres$start.date.date <- strptime(drealpres$start.date, format="%d/%m/%Y")
drealpres$end.date.date <- strptime(drealpres$end.date, format="%d/%m/%Y")
dtemp2 <- data.frame(date=seq(min(drealpres$start.date.date),max(drealpres$end.date.date), by=60*60*24))
dtemp2$dayofyear <- yday(dtemp2$date)
dtemp2 <- dtemp2[dtemp2$dayofyear<=yday(strptime("01.06.2020", format="%d.%m.%Y")),]
for(i in 1:sum(drealpres$animal=="horse")){
ii <- c(1:nrow(drealpres))[drealpres$animal=="horse"][i]
dtemp2$horse[dtemp2$date>=drealpres$start.date.date[ii] & dtemp2$date<=drealpres$end.date.date[ii]] <- drealpres$number[ii]
}
for(i in 1:sum(drealpres$animal=="cattle")){
ii <- c(1:nrow(drealpres))[drealpres$animal=="cattle"][i]
dtemp2$cattle[dtemp2$date>=drealpres$start.date.date[ii] & dtemp2$date<=drealpres$end.date.date[ii]] <- drealpres$number[ii]
}
dtemp2$cattle[is.na(dtemp2$cattle)] <- 0
dtemp2$horse[is.na(dtemp2$horse)] <- 0
dtemp2$year <- year(dtemp2$date)
dtemp3 <- aggregate(cattle~year, data=dtemp2, FUN=sum)
dtemp3$horse <- aggregate(horse~year, data=dtemp2, FUN=sum)$horse
dgrazersumveg$realpresCattle <- dtemp3$cattle[match(dgrazersumveg$year.season, dtemp3$year)]
dgrazersumveg$realpresHorse <- dtemp3$horse[match(dgrazersumveg$year.season, dtemp3$year)]
dgrazersumveg$cattle.corr <- dgrazersumveg$Cattle/dgrazersumveg$nrCattletracked*dgrazersumveg$realpresCattle
dgrazersumveg$horse.corr <-dgrazersumveg$Horse/dgrazersumveg$nrHorsetracked*dgrazersumveg$realpresHorse
dgrazersumveg$grazer.corr <- dgrazersumveg$horse.corr + dgrazersumveg$cattle.corrInformation about the dataset Flora Indicativa (Landolt et al 2010), which is used to assign functional traits to the sampled plant species.
The following variables of the Flora Indicativa are used:
Climate - L (light): 1 deep shade; 2 shade; 3 semi-shade; 4 well lit; 5 full light [additionally we -> sd light because grazers could affect the between-patch variance]
Soil - N (nutrients): 1 very infertile; 2 infertile; 3 medium infertile to medium fertile; 4 fertile; 5 very fertile and over-rich
Life history - MV (mowing /grazing tolerance) 1 bad; 2 little; 3 moderate; 4 good; 5 very good
(the variable name is in two rows in the original excel file)
Information about the dataset of field plant sampling:
The raw dataset contains typos, which will be corrected below (the original txt dataset remains unchanged).
Elymus pungens occurs in the Petite Camargue but no data exist in Flora Helvetica.
Salix x rubens is a hybrid species in the genus Salix. As the hybrid is not present in the Flora Indicativa, we use the attributes of the type species of the genus, Salix alba.
Species that were not determined to the species level -> NA
## cellnr corresponds to plot_nr
dcellnr <- read.table("../data/cells_info_veg-plots-with-forest.txt", header=TRUE, sep="\t")
dveg$cellnr <- dcellnr$cellnr[match(dveg$plot_nr, dcellnr$plot_nr)]
dflora <- read.table("../data/FLORA-INDICATIVA-EN.txt", header=TRUE, skip=1, sep="\t")
dflora$L..light. <- as.numeric(dflora$L..light.) # x turn to NA## Warning: NAs introduced by coercion
dveg$species1 <- trimws(dveg$species1, which="both")
dveg$species2 <- trimws(dveg$species2, which="both")
dveg$species3 <- trimws(dveg$species3, which="both")
# correct species names
index <- dveg$species1=="Calamagrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Calamagrostis epigejos"
index <- dveg$species2=="Calamagrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Calamagrostis epigejos"
index <- dveg$species3=="Calamagrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Calamagrostis epigejos"
index <- dveg$species1=="Calamgrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Calamagrostis epigejos"
index <- dveg$species2=="Calamgrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Calamagrostis epigejos"
index <- dveg$species3=="Calamgrostis epigeios"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Calamagrostis epigejos"
index <- dveg$species1=="Seculigera varia"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Securigera varia"
index <- dveg$species2=="Seculigera varia"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Securigera varia"
index <- dveg$species3=="Seculigera varia"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Securigera varia"
index <- dveg$species1=="Seculigera Varia"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Securigera varia"
index <- dveg$species2=="Seculigera Varia"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Securigera varia"
index <- dveg$species3=="Seculigera Varia"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Securigera varia"
index <- dveg$species1=="Salix rubens"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Salix alba"
index <- dveg$species2=="Salix rubens"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Salix alba"
index <- dveg$species3=="Salix rubens"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Salix alba"
index <- dveg$species1=="Odontites serotinus"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Odontites vernus"
index <- dveg$species2=="Odontites serotinus"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Odontites vernus"
index <- dveg$species3=="Odontites serotinus"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Odontites vernus"
index <- dveg$species1=="Odontites rubra"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Odontites vernus"
index <- dveg$species2=="Odontites rubra"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Odontites vernus"
index <- dveg$species3=="Odontites rubra"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Odontites vernus"
index <- is.element(dveg$species1,
c("Bromus hordaceous", "Bromus hordaeceus", "Bromus hordaceus", "Bromus hordaceous", "Bromes hordaceous")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Bromus hordeaceus"
index <- is.element(dveg$species2,
c("Bromus hordaceous", "Bromus hordaeceus", "Bromus hordaceus", "Bromus hordaceous", "Bromes hordaceous")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Bromus hordeaceus"
index <- is.element(dveg$species3,
c("Bromus hordaceous", "Bromus hordaeceus", "Bromus hordaceus", "Bromus hordaceous", "Bromes hordaceous")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Bromus hordeaceus"
index <- dveg$species1=="Bromes sterilis"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Bromus sterilis"
index <- dveg$species2=="Bromes sterilis"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Bromus sterilis"
index <- dveg$species3=="Bromes sterilis"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Bromus sterilis"
index <- is.element(dveg$species1,
c("Arenaria serpillyfolia", "Arenaria serpillifolia")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Arenaria serpyllifolia"
index <- is.element(dveg$species2,
c("Arenaria serpillyfolia", "Arenaria serpillifolia")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Arenaria serpyllifolia"
index <- is.element(dveg$species3,
c("Arenaria serpillyfolia", "Arenaria serpillifolia")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Arenaria serpyllifolia"
index <- is.element(dveg$species1,
c("Arrhenaterum elatius", "Arrhenatrum elatius")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Arrhenatherum elatius"
index <- is.element(dveg$species2,
c("Arrhenaterum elatius", "Arrhenatrum elatius")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Arrhenatherum elatius"
index <- is.element(dveg$species3,
c("Arrhenaterum elatius", "Arrhenatrum elatius")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Arrhenatherum elatius"
index <- dveg$species1=="Medicago campestre"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Trifolium campestre"
index <- dveg$species2=="Medicago campestre"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Trifolium campestre"
index <- dveg$species3=="Medicago campestre"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Trifolium campestre"
index <- dveg$species1=="Festuca pratensis"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Festuca pratensis"
index <- dveg$species2=="Festuca pratensis"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Festuca pratensis"
index <- dveg$species3=="Festuca pratensis"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Festuca pratensis"
index <- is.element(dveg$species1,
c("APetraraghia prolifera", "Petraraghia prolifera", "Petrarhaghia prolifera", "Petroraghia prolifera")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Petrorhagia prolifera"
index <- is.element(dveg$species2,
c("APetraraghia prolifera", "Petraraghia prolifera", "Petrarhaghia prolifera", "Petroraghia prolifera")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Petrorhagia prolifera"
index <- is.element(dveg$species3,
c("APetraraghia prolifera", "Petraraghia prolifera", "Petrarhaghia prolifera", "Petroraghia prolifera")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Petrorhagia prolifera"
index <- dveg$species1=="Gallium mollugo"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Galium mollugo"
index <- dveg$species2=="Gallium mollugo"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Galium mollugo"
index <- dveg$species3=="Gallium mollugo"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Galium mollugo"
index <- dveg$species1=="Plantage lanceolata"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Plantago lanceolata"
index <- dveg$species2=="Plantage lanceolata"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Plantago lanceolata"
index <- dveg$species3=="Plantage lanceolata"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Plantago lanceolata"
index <- dveg$species1=="Erigeron annus"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Erigeron annuus"
index <- dveg$species2=="Erigeron annus"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Erigeron annuus"
index <- dveg$species3=="Erigeron annus"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Erigeron annuus"
index <- dveg$species1=="Carex sp"; index[is.na(index)] <- FALSE
dveg$species1[index] <- NA
index <- dveg$species2=="Carex sp"; index[is.na(index)] <- FALSE
dveg$species2[index] <- NA
index <- dveg$species3=="Carex sp"; index[is.na(index)] <- FALSE
dveg$species3[index] <- NA
index <- is.element(dveg$species1,
c("Juncus inflxus", "Juncos inflexus")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Juncus inflexus"
index <- is.element(dveg$species2,
c("Juncus inflxus", "Juncos inflexus")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Juncus inflexus"
index <- is.element(dveg$species3,
c("Juncus inflxus", "Juncos inflexus")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Juncus inflexus"
index <- dveg$species1=="Myosotis arvense"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Myosotis arvensis"
index <- dveg$species2=="Myosotis arvense"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Myosotis arvensis"
index <- dveg$species3=="Myosotis arvense"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Myosotis arvensis"
index <- dveg$species1=="Tussilagi farfara"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Tussilago farfara"
index <- dveg$species2=="Tussilagi farfara"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Tussilago farfara"
index <- dveg$species3=="Tussilagi farfara"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Tussilago farfara"
index <- dveg$species1=="Phragmites communis"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Phragmites australis"
index <- dveg$species2=="Phragmites communis"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Phragmites australis"
index <- dveg$species3=="Phragmites communis"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Phragmites australis"
index <- dveg$species1=="Argrostis stolonifera"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Agrostis stolonifera"
index <- dveg$species2=="Argrostis stolonifera"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Agrostis stolonifera"
index <- dveg$species3=="Argrostis stolonifera"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Agrostis stolonifera"
index <- dveg$species1=="Amorpha sp"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Amorpha fruticosa"
index <- dveg$species2=="Amorpha sp"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Amorpha fruticosa"
index <- dveg$species3=="Amorpha sp"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Amorpha fruticosa"
index <- dveg$species1=="Cirsium sp"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Cirsium vulgare"
index <- dveg$species2=="Cirsium sp"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Cirsium vulgare"
index <- dveg$species3=="Cirsium sp"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Cirsium vulgare"
# (this is the most common...)
index <- dveg$species1=="Echyum vulgare"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Echium vulgare"
index <- dveg$species2=="Echyum vulgare"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Echium vulgare"
index <- dveg$species3=="Echyum vulgare"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Echium vulgare"
index <- dveg$species1=="Menta aquatica"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Mentha aquatica"
index <- dveg$species2=="Menta aquatica"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Mentha aquatica"
index <- dveg$species3=="Menta aquatica"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Mentha aquatica"
index <- dveg$species1=="Medication lupulina"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Medicago lupulina"
index <- dveg$species2=="Medication lupulina"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Medicago lupulina"
index <- dveg$species3=="Medication lupulina"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Medicago lupulina"
index <- dveg$species1=="Cinosorus cristatus"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Cynosurus cristatus"
index <- dveg$species2=="Cinosorus cristatus"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Cynosurus cristatus"
index <- dveg$species3=="Cinosorus cristatus"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Cynosurus cristatus"
index <- dveg$species1=="Argostis gigantea"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Agrostis gigantea"
index <- dveg$species2=="Argostis gigantea"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Agrostis gigantea"
index <- dveg$species3=="Argostis gigantea"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Agrostis gigantea"
index <- dveg$species1=="Dactilys glomerata"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Dactylis glomerata"
index <- dveg$species2=="Dactilys glomerata"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Dactylis glomerata"
index <- dveg$species3=="Dactilys glomerata"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Dactylis glomerata"
index <- is.element(dveg$species1,
c("Picris hieraceoides", "Picris hieraceous", "Picris hyeracioides")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Picris hieracioides"
index <- is.element(dveg$species2,
c("Picris hieraceoides", "Picris hieraceous", "Picris hyeracioides")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Picris hieracioides"
index <- is.element(dveg$species3,
c("Picris hieraceoides", "Picris hieraceous", "Picris hyeracioides")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Picris hieracioides"
index <- is.element(dveg$species1,
c("Hypericum officinalis")); index[is.na(index)] <- FALSE
dveg$species1[index] <- "Hypericum perforatum"
index <- is.element(dveg$species2,
c("Hypericum officinalis")); index[is.na(index)] <- FALSE
dveg$species2[index] <- "Hypericum perforatum"
index <- is.element(dveg$species3,
c("Hypericum officinalis")); index[is.na(index)] <- FALSE
dveg$species3[index] <- "Hypericum perforatum"
index <- dveg$species1=="Potentilla reptans ?"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Potentilla reptans"
index <- dveg$species2=="Potentilla reptans ?"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Potentilla reptans"
index <- dveg$species3=="Potentilla reptans ?"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Potentilla reptans"
#dveg[dveg==“Populus balsamifera”] <- NA
# (Flora Indicativa does not contain, and not common in the PCA and in the data)
index <- dveg$species1=="Vicia cracca tenuifolia"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Vicia cracca"
index <- dveg$species2=="Vicia cracca tenuifolia"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Vicia cracca"
index <- dveg$species3=="Vicia cracca tenuifolia"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Vicia cracca"
index <- dveg$species1=="Origanium vulgare"; index[is.na(index)] <- FALSE
dveg$species1[index] <- "Origanum vulgare"
index <- dveg$species2=="Origanium vulgare"; index[is.na(index)] <- FALSE
dveg$species2[index] <- "Origanum vulgare"
index <- dveg$species3=="Origanium vulgare"; index[is.na(index)] <- FALSE
dveg$species3[index] <- "Origanum vulgare"
# end of correction of species namesdveg$Light_sp1 <- dflora$L..light.[match(dveg$species1, dflora$Taxon)]
dveg$Light_sp2 <- dflora$L..light.[match(dveg$species2, dflora$Taxon)]
dveg$Light_sp3 <- dflora$L..light.[match(dveg$species3, dflora$Taxon)]
unique(dveg$species1[is.na(dveg$Light_sp1)])## [1] NA "Elymus pungens" ""
## [1] NA "Elymus pungens" "Ranunculus sp" "Artemisia sp"
## [5] ""
## [1] NA "Populus balsamifera" ""
dveg$N_sp1 <- dflora$N..Nutrients.[match(dveg$species1, dflora$Taxon)]
dveg$N_sp2 <- dflora$N..Nutrients.[match(dveg$species2, dflora$Taxon)]
dveg$N_sp3 <- dflora$N..Nutrients.[match(dveg$species3, dflora$Taxon)]
dveg$MV_sp1 <- dflora$MV..mowing.tolerance.[match(dveg$species1, dflora$Taxon)]
dveg$MV_sp2 <- dflora$MV..mowing.tolerance.[match(dveg$species2, dflora$Taxon)]
dveg$MV_sp3 <- dflora$MV..mowing.tolerance.[match(dveg$species3, dflora$Taxon)]
dveg$light.mean <- apply(dveg[,c("Light_sp1","Light_sp2","Light_sp3")], 1, mean, na.rm=TRUE)
dveg$light.sd <- apply(dveg[,c("Light_sp1","Light_sp2","Light_sp3")], 1, sd, na.rm=TRUE)
dveg$N_sp1 <- as.numeric(dveg$N_sp1)
dveg$N_sp2 <- as.numeric(dveg$N_sp2)
dveg$N_sp3 <- as.numeric(dveg$N_sp3)
dveg$N.mean <- apply(dveg[,c("N_sp1","N_sp2","N_sp3")], 1, mean, na.rm=TRUE)
dveg$MV_sp1 <- as.numeric(dveg$MV_sp1)## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
dveg$MV.mean <- apply(dveg[,c("MV_sp1","MV_sp2","MV_sp3")], 1, mean, na.rm=TRUE)
# add grazerdensity
dveg$horse.corr <- dgrazersumveg$horse.corr[match(paste(dveg$cellnr, dveg$year), paste(dgrazersumveg$cellnr, dgrazersumveg$year.season))]
dveg$cattle.corr <- dgrazersumveg$cattle.corr[match(paste(dveg$cellnr, dveg$year), paste(dgrazersumveg$cellnr, dgrazersumveg$year.season))]
dveg$grazer.corr <- dgrazersumveg$grazer.corr[match(paste(dveg$cellnr, dveg$year), paste(dgrazersumveg$cellnr, dgrazersumveg$year.season))]
dveg$grazer.corr[dveg$year==2018] <- 0
dveg$horse.corr[dveg$year==2018] <- 0
dveg$cattle.corr[dveg$year==2018] <- 0
dveg <- dveg[!is.na(dveg$year),]
dveg$heightm <- apply(dveg[,c("height1_1", "height1_2", "height1_3" ,"height1_4" , "height1_5" , "height1_6", "height1_7" , "height1_8" , "height1_9" ,
"height1_10" , "height2_1","height2_2", "height2_3", "height2_4", "height2_5",
"height2_6", "height2_7", "height2_8", "height2_9", "height2_10")], 1, mean)
pairs(dveg[,c("grazer.corr", "horse.corr", "cattle.corr", "light.mean", "light.sd", "N.mean", "MV.mean", "heightm", "total_cover")], panel=panel.smooth)mod <- lmer(grazer.log~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)) + (1|cellnr), data=dat)## Warning: Some predictor variables are on very different scales: consider
## rescaling
resdat <- data.frame(resid=resid(mod), x=dcells$x[match(dat$cellnr, dcells$cellnr)],
y=dcells$y[match(dat$cellnr, dcells$cellnr)])
variomod <- vgram(resdat[,c("x", "y")], resdat$resid,
type="variogram")
par(mfrow=c(1,2))
plot(variomod, ylim=c(0,0.8))
# bubble plot
plot(jitter(dat$x,20), jitter(dat$y, 20), pch=16, col=c("blue", "orange")[(sign(resid(mod))+3)/2], cex=abs(resid(mod)))Figure 3.1: Semi-variogram and bubble plot of the residuals of the grazer model to check for spatial correlation. Positive residuals are orange, negative residuals are blue.
modh <- lmer(horse.log~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)) + (1|cellnr), data=dat)## Warning: Some predictor variables are on very different scales: consider
## rescaling
resdat <- data.frame(resid=resid(modh), x=dcells$x[match(dat$cellnr, dcells$cellnr)],
y=dcells$y[match(dat$cellnr, dcells$cellnr)])
variomod <- vgram(resdat[,c("x", "y")], resdat$resid,
type="variogram")
par(mfrow=c(1,2))
plot(variomod, ylim=c(0,0.8))
# bubble plot
plot(jitter(dat$x,20), jitter(dat$y, 20), pch=16, col=c("blue", "orange")[(sign(resid(modh))+3)/2], cex=abs(resid(mod)))Figure 3.2: Visualisation of the residuals of the horse model to check for spatial correlation. Positive residuals are orange, negative residuals are blue.
modc <- lmer(cattle.log~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)) + (1|cellnr), data=dat)## Warning: Some predictor variables are on very different scales: consider
## rescaling
resdat <- data.frame(resid=resid(modc), x=dcells$x[match(dat$cellnr, dcells$cellnr)],
y=dcells$y[match(dat$cellnr, dcells$cellnr)])
variomod <- vgram(resdat[,c("x", "y")], resdat$resid,
type="variogram")
par(mfrow=c(1,2))
plot(variomod)
# bubble plot
plot(jitter(dat$x,20), jitter(dat$y, 20), pch=16, col=c("blue", "orange")[(sign(resid(modc))+3)/2], cex=abs(resid(mod)))Figure 3.3: Visualisation of the residuals of the cattle model to check for spatial correlation. Positive residuals are orange, negative residuals are blue.
No obvious spatial correlation can be seen in the residuals, i.e., both negative and positive residuals are spread around the whole study area. However, in the eastern part of the study area, the absolute values of the residuals are larger. This pattern indicates that in the eastern part of the study area the variance in grazing intensity is larger compared to the western part. When plotting the residuals against the fitted values (plot not shown), we see that the large absolute residuals are for low fitted values. Therefore, this may be a consequence of low sample size and the logarithm transformation. I would accept that heterogeneity of variance.
dvegc <- dveg[complete.cases(dveg[, c("light.mean", "light.sd", "N.mean", "MV.mean", "horse.corr", "cattle.corr", "year", "cellnr", "heightm", "total_cover")]),]
dvegc$year.z <- as.numeric(scale(dvegc$year))
cor(dvegc[, c("light.mean", "light.sd", "N.mean", "MV.mean", "heightm", "total_cover")])## light.mean light.sd N.mean MV.mean heightm
## light.mean 1.00000000 0.21296431 0.163526150 0.032236377 -0.15073519
## light.sd 0.21296431 1.00000000 0.113114205 -0.040651977 0.01498528
## N.mean 0.16352615 0.11311420 1.000000000 -0.008692668 0.19598873
## MV.mean 0.03223638 -0.04065198 -0.008692668 1.000000000 -0.33577522
## heightm -0.15073519 0.01498528 0.195988730 -0.335775224 1.00000000
## total_cover -0.06599680 0.12281765 0.020094808 0.029852613 0.33171901
## total_cover
## light.mean -0.06599680
## light.sd 0.12281765
## N.mean 0.02009481
## MV.mean 0.02985261
## heightm 0.33171901
## total_cover 1.00000000
## [1] 11.15824
## [1] 5.080214
dvegc$horse.corr.sqrt <- sqrt(dvegc$horse.corr)
dvegc$cattle.corr.sqrt <- sqrt(dvegc$cattle.corr)
modlightm <- lmer(light.mean ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modlightsd <- lmer(light.sd ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modN.m <- lmer(N.mean ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modMV.m <- lmer(MV.mean ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modheightm <- lmer(log(heightm) ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
modcover <- lmer(total_cover ~ horse.corr.sqrt + cattle.corr.sqrt + year.z + (1|cellnr), data=dvegc)
# plot(modlightm)
# plot(modlightsd)
# plot(modN.m)
# plot(modMV.m)
# plot(modheightm)
# plot(modcover)
# because data is ordered according to cellnr within year,
# the following plots should show spatial correlation if present
# acf(resid(modlightm))
# acf(resid(modlightsd))
# acf(resid(modN.m))
# acf(resid(modMV.m))
# acf(resid(modheightm))
# acf(resid(modcover))
# -> no spatial correlation
dvegc$x <- dcells$x[match(dvegc$cellnr, dcells$cellnr)]
dvegc$y <- dcells$y[match(dvegc$cellnr, dcells$cellnr)]
resdat <- data.frame(residLm=resid(modlightm),
residLsd=resid(modlightsd),
residNm=resid(modN.m),
residMVm=resid(modMV.m),
residheight=resid(modheightm),
residcover=resid(modcover),
x=dvegc$x,
y=dvegc$y)
varLM <- vgram(resdat[,c("x", "y")], resdat$residLm,
type="variogram")
varLsd <- vgram(resdat[,c("x", "y")], resdat$residLsd,
type="variogram")
varNm <- vgram(resdat[,c("x", "y")], resdat$residNm,
type="variogram")
varMVm <- vgram(resdat[,c("x", "y")], resdat$residMVm,
type="variogram")
varheightm <- vgram(resdat[,c("x", "y")], resdat$residheight,
type="variogram")
varcover <- vgram(resdat[,c("x", "y")], resdat$residcover,
type="variogram")
# par(mfrow=c(3,2))
# plot(varLM)
# plot(varLsd)
# plot(varNm)
# plot(varMVm)
# plot(varheightm)
# plot(varcover)
# -> again, no spatial correlation
nsim <- 5000
bsimLm <- sim(modlightm, n.sim=nsim)
bsimLsd <- sim(modlightsd, n.sim=nsim)
bsimNm <- sim(modN.m, n.sim=nsim)
bsimMVm <- sim(modMV.m, n.sim=nsim)
bsimheightm <- sim(modheightm, n.sim=nsim)
bsimcover <- sim(modcover, n.sim=nsim)## [1] "Plot.nr" "cellnr" "Cover.tree" "Cover.sapling"
## [5] "Cover.shrub" "Cover.meadow" "Surface.water" "Bareground"
## 'data.frame': 163 obs. of 8 variables:
## $ Plot.nr : chr "A75" "" "A73" "A74" ...
## $ cellnr : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Cover.tree : int 0 40 0 0 0 30 0 0 0 0 ...
## $ Cover.sapling: int 3 0 1 2 15 5 1 20 12 1 ...
## $ Cover.shrub : int 0 5 0 1 0 0 2 3 3 0 ...
## $ Cover.meadow : int 97 40 99 97 60 60 97 77 35 69 ...
## $ Surface.water: int 0 0 0 0 20 5 0 0 50 30 ...
## $ Bareground : int 0 15 0 0 5 30 0 0 0 0 ...
landcover_cols <- c("Cover.tree", "Cover.sapling", "Cover.shrub", "Cover.meadow", "Surface.water", "Bareground")
landcover_data <- dlandcover[, landcover_cols]
landcover_data <- as.data.frame(lapply(landcover_data, function(x) as.numeric(as.character(x))))
# convert landcover_data for ggplot2
landcover_long <- melt(landcover_data, variable.name = "LandCoverType", value.name = "Percentage")## No id variables; using all as measure variables
# mean and sd for each land cover type
landcover_means <- aggregate(Percentage ~ LandCoverType, data = landcover_long, FUN = mean)
landcover_sds <- aggregate(Percentage ~ LandCoverType, data = landcover_long, FUN = sd)
landcover_rename <- c(
"Cover.tree" = "tree",
"Cover.sapling" = "sapling",
"Cover.shrub" = "shrub",
"Cover.meadow" = "meadow",
"Surface.water" = "water surface",
"Bareground" = "bareground"
)
landcover_long$LandCoverType <- landcover_rename[as.character(landcover_long$LandCoverType)]
landcover_means$LandCoverType <- landcover_rename[as.character(landcover_means$LandCoverType)]
landcover_sds$LandCoverType <- landcover_rename[as.character(landcover_sds$LandCoverType)]
landcover_colors <- c(
"bareground" = "#D6BD8DFF",
"meadow" = "#F6EDBDFF",
"sapling" = "#B5B991FF",
"shrub" = "#778868FF",
"tree" = "#3D5941FF",
"water surface" = "#2887A1FF"
)
landcover_long$LandCoverType <- factor(landcover_long$LandCoverType, levels = names(landcover_colors))
# plot
p <- ggplot(landcover_long, aes(x = LandCoverType, y = Percentage, fill = LandCoverType)) +
geom_violin(trim = TRUE, scale = "width", adjust = 1.2) + # adjust the scale and width of the violins
geom_point(data = landcover_means, aes(x = LandCoverType, y = Percentage), color = "white", size = 3) +
geom_errorbar(data = landcover_means,
aes(x = LandCoverType,
ymin = pmax(0, Percentage - landcover_sds$Percentage),
ymax = Percentage + landcover_sds$Percentage),
color = "white", width = 0.1) +
geom_text(data = landcover_means, aes(x = LandCoverType, y = Percentage, label = round(Percentage, 1)), color = "white", vjust = 1, hjust = -0.4, size = 3.5) +
scale_fill_manual(values = landcover_colors) +
labs(title = "Macrohabitat composition on the Rhine Island", x = "Land cover type", y = "Percentage") +
theme_minimal(base_size = 15) +
theme(legend.position = "none", panel.background = element_rect(fill = "lightgrey"))
pnsim <- 5000
bsim <- sim(mod, n.sim=nsim)
tab <- t(apply(bsim@fixef, 2, quantile, prob=c(0.5, 0.025, 0.975)))
nfixef <- nrow(tab)
# add random effects
sdcellnr <- apply(bsim@ranef$cellnr[,,1], 1, sd)
tab <- rbind(tab, quantile(sdcellnr, probs=c(0.5, 0.025, 0.975)))
tab <- rbind(tab, quantile(bsim@sigma, probs=c(0.5, 0.025, 0.975)))
rownames(tab)[c(nfixef+1, nfixef+2)] <- c("cellnr_sd", "sigma")
kable(tab, dig=2, caption="Median, 2.5% and 97.5% quantiles of the posterior distributions of the model coefficients for the model with grazer density as response.")| 50% | 2.5% | 97.5% | |
|---|---|---|---|
| (Intercept) | -1.81 | -2.21 | -1.43 |
| Cover.tree | 0.00 | -0.02 | 0.02 |
| Cover.sapling | 0.00 | -0.02 | 0.03 |
| Cover.shrub | -0.09 | -0.19 | 0.02 |
| Surface.water | 0.04 | 0.01 | 0.08 |
| Bareground | 0.02 | -0.03 | 0.07 |
| I(Cover.tree^2) | 0.00 | 0.00 | 0.00 |
| I(Cover.sapling^2) | 0.00 | 0.00 | 0.00 |
| I(Cover.shrub^2) | 0.01 | 0.00 | 0.02 |
| I(Surface.water^2) | 0.00 | 0.00 | 0.00 |
| I(Bareground^2) | 0.00 | 0.00 | 0.00 |
| seasonwinter | -0.27 | -0.53 | -0.02 |
| Cover.tree:seasonwinter | -0.02 | -0.03 | 0.00 |
| Cover.sapling:seasonwinter | 0.02 | 0.01 | 0.04 |
| Cover.shrub:seasonwinter | 0.18 | 0.11 | 0.25 |
| Surface.water:seasonwinter | -0.03 | -0.05 | -0.01 |
| Bareground:seasonwinter | 0.00 | -0.04 | 0.03 |
| I(Cover.tree^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| I(Cover.sapling^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| I(Cover.shrub^2):seasonwinter | -0.02 | -0.02 | -0.01 |
| I(Surface.water^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| I(Bareground^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| cellnr_sd | 0.83 | 0.78 | 0.88 |
| sigma | 0.76 | 0.72 | 0.79 |
mprop <- apply(dat[,cov.varnames], 2, mean)
par(mfrow=c(3,2), mar=c(4,3,0.5,0.5), mgp=c(2,0.5,0), oma=c(0, 0, 4, 0))
for(cov.var in cov.varnames){
range.var <- range(dat[,cov.var])
eval(parse(text=paste0("newdat <- expand.grid(",cov.var,"=seq(range.var[1], range.var[2], length=100),
season=levels(factor(dat$season)))")))
restvar <- cov.varnames[!is.element(cov.varnames, cov.var)]
for(j in restvar){
newdat[[j]] <- mprop[j]/(100-mprop[cov.var])*(100-newdat[,cov.var])
}
#Xmat <- model.matrix(~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + season + season:(Cover.tree + Cover.sapling + Cover.shrub + Surface.water+ Bareground), data=newdat)
Xmat <- model.matrix(~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground +
I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)), data=newdat)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- Xmat%*%bsim@fixef[i,]
newdat$fit <- apply(fitmat, 1, mean)
newdat$lwr <- apply(fitmat, 1, quantile, probs=0.025)
newdat$upr <- apply(fitmat, 1, quantile, probs=0.975)
ylimforlater <- range(dat$grazer.log)
plot(jitter(dat[,cov.var]), dat$grazer.log, xlab="", ylab="log(sum of locations)", col=c("#FFC20A", "#0C7BDC")[as.numeric(dat$season)])
title(xlab=paste(cov.var, "[%]"), cex.lab=1.3)
index <- newdat$season=="summer"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#FFC20A")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#FFC20A")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#FFC20A")
index <- newdat$season=="winter"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#0C7BDC")
mtext("Grazer density", side=3, line=1, outer=TRUE, cex=1.5)
}# close cov.varFigure 4.1: Grazer density in relation to the cover variables for summer (orange) and winter (blue). 95% uncertainty intervals are given.
nsim <- 5000
bsim <- sim(modh, n.sim=nsim)
tab <- t(apply(bsim@fixef, 2, quantile, prob=c(0.5, 0.025, 0.975)))
nfixef <- nrow(tab)
# add random effects
sdcellnr <- apply(bsim@ranef$cellnr[,,1], 1, sd)
tab <- rbind(tab, quantile(sdcellnr, probs=c(0.5, 0.025, 0.975)))
tab <- rbind(tab, quantile(bsim@sigma, probs=c(0.5, 0.025, 0.975)))
rownames(tab)[c(nfixef+1, nfixef+2)] <- c("cellnr_sd", "sigma")
kable(tab, dig=2, caption="Median, 2.5% and 97.5% quantiles of the posterior distributions of the model coefficients for the model with horse location density as response.")| 50% | 2.5% | 97.5% | |
|---|---|---|---|
| (Intercept) | -2.46 | -2.84 | -2.06 |
| Cover.tree | 0.00 | -0.02 | 0.02 |
| Cover.sapling | 0.00 | -0.02 | 0.03 |
| Cover.shrub | -0.11 | -0.22 | -0.01 |
| Surface.water | 0.02 | -0.01 | 0.05 |
| Bareground | 0.02 | -0.03 | 0.07 |
| I(Cover.tree^2) | 0.00 | 0.00 | 0.00 |
| I(Cover.sapling^2) | 0.00 | 0.00 | 0.00 |
| I(Cover.shrub^2) | 0.01 | 0.00 | 0.02 |
| I(Surface.water^2) | 0.00 | 0.00 | 0.00 |
| I(Bareground^2) | 0.00 | 0.00 | 0.00 |
| seasonwinter | -0.13 | -0.37 | 0.11 |
| Cover.tree:seasonwinter | -0.01 | -0.02 | 0.00 |
| Cover.sapling:seasonwinter | 0.03 | 0.01 | 0.04 |
| Cover.shrub:seasonwinter | 0.19 | 0.13 | 0.26 |
| Surface.water:seasonwinter | 0.00 | -0.02 | 0.02 |
| Bareground:seasonwinter | -0.01 | -0.04 | 0.03 |
| I(Cover.tree^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| I(Cover.sapling^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| I(Cover.shrub^2):seasonwinter | -0.02 | -0.02 | -0.01 |
| I(Surface.water^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| I(Bareground^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| cellnr_sd | 0.83 | 0.79 | 0.88 |
| sigma | 0.69 | 0.66 | 0.72 |
par(mfrow=c(3,2), mar=c(4,3,0.5,0.5), mgp=c(2,0.5,0), oma=c(0, 0, 4, 0))
for(cov.var in cov.varnames){
range.var <- range(dat[,cov.var])
eval(parse(text=paste0("newdat <- expand.grid(",cov.var,"=seq(range.var[1], range.var[2], length=100),
season=levels(factor(dat$season)))")))
restvar <- cov.varnames[!is.element(cov.varnames, cov.var)]
for(j in restvar){
newdat[[j]] <- mprop[j]/(100-mprop[cov.var])*(100-newdat[,cov.var])
}
Xmat <- model.matrix(~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground +
I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)), data=newdat)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- Xmat%*%bsim@fixef[i,]
newdat$fit <- apply(fitmat, 1, mean)
newdat$lwr <- apply(fitmat, 1, quantile, probs=0.025)
newdat$upr <- apply(fitmat, 1, quantile, probs=0.975)
plot(jitter(dat[,cov.var]), dat$horse.log, xlab="", ylab="log(sum of locations)", col=c("#FFC20A", "#0C7BDC")[as.numeric(dat$season)], ylim=ylimforlater)
title(xlab=paste(cov.var, "[%]"), cex.lab=1.3)
index <- newdat$season=="summer"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#FFC20A")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#FFC20A")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#FFC20A")
index <- newdat$season=="winter"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#0C7BDC")
mtext("Horse density", side=3, line=1, outer=TRUE, cex=1.5)
}# close cov.varFigure 4.2: Horse location density in relation to the cover variables for summer (orange) and winter (blue). 95% uncertainty intervals are given.
nsim <- 5000
bsim <- sim(modc, n.sim=nsim)
tab <- t(apply(bsim@fixef, 2, quantile, prob=c(0.5, 0.025, 0.975)))
nfixef <- nrow(tab)
# add random effects
sdcellnr <- apply(bsim@ranef$cellnr[,,1], 1, sd)
tab <- rbind(tab, quantile(sdcellnr, probs=c(0.5, 0.025, 0.975)))
tab <- rbind(tab, quantile(bsim@sigma, probs=c(0.5, 0.025, 0.975)))
rownames(tab)[c(nfixef+1, nfixef+2)] <- c("cellnr_sd", "sigma")
kable(tab, dig=2, caption="Median, 2.5% and 97.5% quantiles of the posterior distributions of the model coefficients for the model with horse location density as response.")| 50% | 2.5% | 97.5% | |
|---|---|---|---|
| (Intercept) | -2.63 | -2.98 | -2.27 |
| Cover.tree | 0.00 | -0.02 | 0.02 |
| Cover.sapling | 0.00 | -0.02 | 0.03 |
| Cover.shrub | -0.04 | -0.13 | 0.06 |
| Surface.water | 0.06 | 0.03 | 0.09 |
| Bareground | 0.01 | -0.03 | 0.06 |
| I(Cover.tree^2) | 0.00 | 0.00 | 0.00 |
| I(Cover.sapling^2) | 0.00 | 0.00 | 0.00 |
| I(Cover.shrub^2) | 0.01 | 0.00 | 0.01 |
| I(Surface.water^2) | 0.00 | 0.00 | 0.00 |
| I(Bareground^2) | 0.00 | 0.00 | 0.00 |
| seasonwinter | -0.82 | -1.18 | -0.47 |
| Cover.tree:seasonwinter | -0.01 | -0.03 | 0.01 |
| Cover.sapling:seasonwinter | 0.02 | 0.00 | 0.05 |
| Cover.shrub:seasonwinter | 0.10 | 0.00 | 0.19 |
| Surface.water:seasonwinter | -0.05 | -0.08 | -0.02 |
| Bareground:seasonwinter | 0.01 | -0.04 | 0.06 |
| I(Cover.tree^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| I(Cover.sapling^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| I(Cover.shrub^2):seasonwinter | -0.01 | -0.02 | 0.00 |
| I(Surface.water^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| I(Bareground^2):seasonwinter | 0.00 | 0.00 | 0.00 |
| cellnr_sd | 0.59 | 0.54 | 0.64 |
| sigma | 1.04 | 1.00 | 1.09 |
par(mfrow=c(3,2), mar=c(4,3,0.5,0.5), mgp=c(2,0.5,0), oma=c(0, 0, 4, 0))
for(cov.var in cov.varnames){
range.var <- range(dat[,cov.var])
eval(parse(text=paste0("newdat <- expand.grid(",cov.var,"=seq(range.var[1], range.var[2], length=100),
season=levels(factor(dat$season)))")))
restvar <- cov.varnames[!is.element(cov.varnames, cov.var)]
for(j in restvar){
newdat[[j]] <- mprop[j]/(100-mprop[cov.var])*(100-newdat[,cov.var])
}
Xmat <- model.matrix(~Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2) + season + season:(Cover.tree + Cover.sapling + Cover.shrub+ Surface.water+ Bareground + I(Cover.tree^2) + I(Cover.sapling^2) + I(Cover.shrub^2)+ I(Surface.water^2)+ I(Bareground^2)), data=newdat)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
for(i in 1:nsim) fitmat[,i] <- Xmat%*%bsim@fixef[i,]
newdat$fit <- apply(fitmat, 1, mean)
newdat$lwr <- apply(fitmat, 1, quantile, probs=0.025)
newdat$upr <- apply(fitmat, 1, quantile, probs=0.975)
plot(jitter(dat[,cov.var]), dat$cattle.log, xlab="", ylab="log(sum of locations)", col=c("#FFC20A", "#0C7BDC")[as.numeric(dat$season)], ylim=ylimforlater)
title(xlab=paste(cov.var, "[%]"), cex.lab=1.3)
index <- newdat$season=="summer"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#FFC20A")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#FFC20A")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#FFC20A")
index <- newdat$season=="winter"
lines(newdat[index, cov.var], newdat$fit[index], lwd=2, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$lwr[index], lwd=1, lty=3, col="#0C7BDC")
lines(newdat[index, cov.var], newdat$upr[index], lwd=1, lty=3, col="#0C7BDC")
mtext("Cattle density", side=3, line=1, outer=TRUE, cex=1.5)
}# close cov.varFigure 4.3: Cattle location density in relation to the cover variables for summer (orange) and winter (blue). 95% uncertainty intervals are given.
Diversity of vegetation on the study site in the different years. The diversity indices are measured taking the number of occupied cells per species as if they were number of individuals.
dveglong <- data.frame(species=c(dveg$species1, dveg$species2, dveg$species3),
year=rep(dveg$year, 3),
season=rep(dveg$season,3),
plot_nr=rep(dveg$plot_nr,3))
ddiv <- aggregate(plot_nr~species + year+ season, data=dveglong, FUN=length)
ddiversity <- data.frame(year=2018:2021,
nrspecies = as.numeric(table(ddiv$year)),
totalcells=tapply(ddiv$plot_nr, ddiv$year, FUN=sum))
ddiv$totalcells <- ddiversity$totalcells[match(ddiv$year, ddiversity$year)]
ddiv$p <- ddiv$plot_nr/ddiv$totalcells
shannon <- function(p) -sum(p*log(p)) # log is natural logarithm
ddiversity$Shannon <- tapply(ddiv$p, ddiv$year, FUN=shannon)
ddiversity$Pielou_evenness <- ddiversity$Shannon/log(ddiversity$nrspecies)
kable(ddiversity[,-3], dig=3, caption="Species richness (nrspecies), Shannon diversity index and Pielou evenness for each year") | year | nrspecies | Shannon | Pielou_evenness | |
|---|---|---|---|---|
| 2018 | 2018 | 51 | 3.352321 | 0.8526119 |
| 2019 | 2019 | 46 | 3.263625 | 0.8524238 |
| 2020 | 2020 | 42 | 3.099439 | 0.8292438 |
| 2021 | 2021 | 43 | 3.138654 | 0.8344821 |
| Interc. | Int.lwr | Int.upr | horse | horse.lwr | horse.upr | cattle | cattle.lwr | cattle.upr | year.z | year.lwr | year.upr | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Light.mean | 3.38 | 3.30 | 3.46 | 0.00 | 0.00 | 0.01 | 0.00 | 0.00 | 0.01 | -0.03 | -0.07 | 0.01 |
| Light.SD | 0.40 | 0.33 | 0.47 | 0.00 | -0.01 | 0.01 | 0.01 | 0.00 | 0.02 | -0.02 | -0.06 | 0.01 |
| N.mean | 3.02 | 2.90 | 3.15 | 0.00 | -0.02 | 0.01 | 0.01 | -0.01 | 0.02 | -0.07 | -0.13 | -0.01 |
| MV.mean | 2.89 | 2.73 | 3.04 | 0.01 | 0.00 | 0.03 | -0.02 | -0.04 | 0.00 | 0.23 | 0.15 | 0.31 |
| log(height.m) | 2.45 | 2.28 | 2.62 | -0.01 | -0.02 | 0.01 | -0.01 | -0.03 | 0.01 | -0.25 | -0.34 | -0.17 |
| total cover | 70.65 | 64.63 | 76.72 | 0.67 | 0.15 | 1.22 | -0.30 | -0.85 | 0.24 | -0.06 | -2.55 | 2.44 |
newdathorse <- data.frame(horse.corr=seq(0, 510, length=100),
cattle.corr=mean(dvegc$cattle.corr),
year.z=0)
newdathorse$horse.corr.sqrt <- sqrt(newdathorse$horse.corr)
newdathorse$cattle.corr.sqrt <- sqrt(newdathorse$cattle.corr)
newdatcattle <- data.frame(cattle.corr=seq(0, 800, length=100),
horse.corr=mean(dvegc$horse.corr),
year.z=0)
newdatcattle$horse.corr.sqrt <- sqrt(newdatcattle$horse.corr)
newdatcattle$cattle.corr.sqrt <- sqrt(newdatcattle$cattle.corr)
newdatyear <- data.frame(cattle.corr=mean(dvegc$cattle.corr),
horse.corr=mean(dvegc$horse.corr),
year=2018:2021)
newdatyear$horse.corr.sqrt <- sqrt(newdatyear$horse.corr)
newdatyear$cattle.corr.sqrt <- sqrt(newdatyear$cattle.corr)
newdatyear$year.z <- (newdatyear$year-mean(dvegc$year))/sd(dvegc$year)
Xmathorse <- model.matrix(~horse.corr.sqrt + cattle.corr.sqrt + year.z, data=newdathorse)
Xmatcattle <- model.matrix(~horse.corr.sqrt + cattle.corr.sqrt + year.z, data=newdatcattle)
Xmatyear <- model.matrix(~horse.corr.sqrt + cattle.corr.sqrt + year.z, data=newdatyear)
# Light mean
newdathorse$Lm <- predict(modlightm, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimLm@fixef[i,]
newdathorse$Lm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$Lm.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatcattle$Lm <- predict(modlightm, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimLm@fixef[i,]
newdatcattle$Lm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$Lm.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatyear$Lm <- predict(modlightm, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimLm@fixef[i,]
newdatyear$Lm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$Lm.upr <- apply(fitmat, 1, quantile, prob=0.975)
# Light SD
newdathorse$Lsd <- predict(modlightsd, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimLsd@fixef[i,]
newdathorse$Lsd.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$Lsd.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatcattle$Lsd <- predict(modlightsd, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimLsd@fixef[i,]
newdatcattle$Lsd.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$Lsd.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatyear$Lsd <- predict(modlightsd, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimLsd@fixef[i,]
newdatyear$Lsd.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$Lsd.upr <- apply(fitmat, 1, quantile, prob=0.975)
# N mean
newdathorse$Nm <- predict(modN.m, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimNm@fixef[i,]
newdathorse$Nm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$Nm.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatcattle$Nm <- predict(modN.m, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimNm@fixef[i,]
newdatcattle$Nm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$Nm.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatyear$Nm <- predict(modN.m, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimNm@fixef[i,]
newdatyear$Nm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$Nm.upr <- apply(fitmat, 1, quantile, prob=0.975)
#MV mean
newdathorse$MVm <- predict(modMV.m, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimMVm@fixef[i,]
newdathorse$MVm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$MVm.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatcattle$MVm <- predict(modMV.m, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimMVm@fixef[i,]
newdatcattle$MVm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$MVm.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatyear$MVm <- predict(modMV.m, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimMVm@fixef[i,]
newdatyear$MVm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$MVm.upr <- apply(fitmat, 1, quantile, prob=0.975)
#height mean
newdathorse$heightm <- predict(modheightm, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimheightm@fixef[i,]
newdathorse$heightm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$heightm.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatcattle$heightm <- predict(modheightm, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimheightm@fixef[i,]
newdatcattle$heightm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$heightm.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatyear$heightm <- predict(modheightm, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimheightm@fixef[i,]
newdatyear$heightm.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$heightm.upr <- apply(fitmat, 1, quantile, prob=0.975)
#cover mean (%)
newdathorse$cover <- predict(modcover, newdata=newdathorse, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdathorse))
for(i in 1:nsim) fitmat[,i] <- Xmathorse%*%bsimcover@fixef[i,]
newdathorse$cover.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdathorse$cover.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatcattle$cover <- predict(modcover, newdata=newdatcattle, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatcattle))
for(i in 1:nsim) fitmat[,i] <- Xmatcattle%*%bsimcover@fixef[i,]
newdatcattle$cover.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatcattle$cover.upr <- apply(fitmat, 1, quantile, prob=0.975)
newdatyear$cover <- predict(modcover, newdata=newdatyear, re.form=NA)
fitmat <- matrix(ncol=nsim, nrow=nrow(newdatyear))
for(i in 1:nsim) fitmat[,i] <- Xmatyear%*%bsimcover@fixef[i,]
newdatyear$cover.lwr <- apply(fitmat, 1, quantile, prob=0.025)
newdatyear$cover.upr <- apply(fitmat, 1, quantile, prob=0.975)
#light mean
par(mfrow=c(6,3), mar=c(2,2,0.5,0.5), oma=c(0,2,2,0))
plot(dvegc$horse.corr, dvegc$light.mean, pch=16, col="#2C6B4F66", cex=0.5, main=NA, ylab="Light mean", xpd=NA, xlab=NA)
mtext("Horse density", side=3, line=0.5)
lines(newdathorse$horse.corr, newdathorse$Lm, lwd=2)
lines(newdathorse$horse.corr, newdathorse$Lm.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$Lm.upr, lwd=1,lty=3)
plot(dvegc$cattle.corr, dvegc$light.mean, pch=16, col="#2C6B4F66", cex=0.5, main=NA)
mtext("Cattle density", side=3, line=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$Lm, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$Lm.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$Lm.upr, lwd=1,lty=3)
plot(jitter(dvegc$year), dvegc$light.mean, pch=16, col="#2C6B4F66", cex=0.5, main=NA)
mtext("Year", side=3, line=0.5)
lines(newdatyear$year, newdatyear$Lm, lwd=2)
lines(newdatyear$year, newdatyear$Lm.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$Lm.upr, lwd=1,lty=3)
# light sd
plot(dvegc$horse.corr, dvegc$light.sd, pch=16, col="#2C6B4F66", cex=0.5, ylab="Light SD", xpd=NA, xlab="")
lines(newdathorse$horse.corr, newdathorse$Lsd, lwd=2)
lines(newdathorse$horse.corr, newdathorse$Lsd.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$Lsd.upr, lwd=1,lty=3)
plot(dvegc$cattle.corr, dvegc$light.sd, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$Lsd, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$Lsd.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$Lsd.upr, lwd=1,lty=3)
plot(jitter(dvegc$year), dvegc$light.sd, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, newdatyear$Lsd, lwd=2)
lines(newdatyear$year, newdatyear$Lsd.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$Lsd.upr, lwd=1,lty=3)
# N mean
plot(dvegc$horse.corr, dvegc$N.mean, pch=16, col="#2C6B4F66", cex=0.5, ylab="N mean", xpd=NA, xlab=NA)
lines(newdathorse$horse.corr, newdathorse$Nm, lwd=2)
lines(newdathorse$horse.corr, newdathorse$Nm.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$Nm.upr, lwd=1,lty=3)
plot(dvegc$cattle.corr, dvegc$N.mean, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$Nm, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$Nm.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$Nm.upr, lwd=1,lty=3)
plot(jitter(dvegc$year), dvegc$N.mean, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, newdatyear$Nm, lwd=2)
lines(newdatyear$year, newdatyear$Nm.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$Nm.upr, lwd=1,lty=3)
# MV mean
plot(dvegc$horse.corr, dvegc$MV.mean, pch=16, col="#2C6B4F66", cex=0.5, ylab="MV mean", xpd=NA, xlab="")
lines(newdathorse$horse.corr, newdathorse$MVm, lwd=2)
lines(newdathorse$horse.corr, newdathorse$MVm.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$MVm.upr, lwd=1,lty=3)
plot(dvegc$cattle.corr, dvegc$MV.mean, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$MVm, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$MVm.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$MVm.upr, lwd=1,lty=3)
plot(jitter(dvegc$year), dvegc$MV.mean, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, newdatyear$MVm, lwd=2)
lines(newdatyear$year, newdatyear$MVm.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$MVm.upr, lwd=1,lty=3)
# height m
plot(dvegc$horse.corr, dvegc$heightm, pch=16, col="#2C6B4F66", cex=0.5, ylab="Height mean", xpd=NA, xlab="")
lines(newdathorse$horse.corr, exp(newdathorse$heightm), lwd=2)
lines(newdathorse$horse.corr, exp(newdathorse$heightm.lwr), lwd=1,lty=3)
lines(newdathorse$horse.corr, exp(newdathorse$heightm.upr), lwd=1,lty=3)
plot(dvegc$cattle.corr, dvegc$heightm, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, exp(newdatcattle$heightm), lwd=2)
lines(newdatcattle$cattle.corr, exp(newdatcattle$heightm.lwr), lwd=1,lty=3)
lines(newdatcattle$cattle.corr, exp(newdatcattle$heightm.upr), lwd=1,lty=3)
plot(jitter(dvegc$year), dvegc$heightm, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, exp(newdatyear$heightm), lwd=2)
lines(newdatyear$year, exp(newdatyear$heightm.lwr), lwd=1,lty=3)
lines(newdatyear$year, exp(newdatyear$heightm.upr), lwd=1,lty=3)
# cover
plot(dvegc$horse.corr, dvegc$total_cover, pch=16, col="#2C6B4F66", cex=0.5, ylab="Total cover", xpd=NA, xlab="")
lines(newdathorse$horse.corr, newdathorse$cover, lwd=2)
lines(newdathorse$horse.corr, newdathorse$cover.lwr, lwd=1,lty=3)
lines(newdathorse$horse.corr, newdathorse$cover.upr, lwd=1,lty=3)
plot(dvegc$cattle.corr, dvegc$total_cover, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatcattle$cattle.corr, newdatcattle$cover, lwd=2)
lines(newdatcattle$cattle.corr, newdatcattle$cover.lwr, lwd=1,lty=3)
lines(newdatcattle$cattle.corr, newdatcattle$cover.upr, lwd=1,lty=3)
plot(jitter(dvegc$year), dvegc$total_cover, pch=16, col="#2C6B4F66", cex=0.5)
lines(newdatyear$year, newdatyear$cover, lwd=2)
lines(newdatyear$year, newdatyear$cover.lwr, lwd=1,lty=3)
lines(newdatyear$year, newdatyear$cover.upr, lwd=1,lty=3)Figure 4.4: Relationship between horse and cattle densities, and plant functional and structural traits